nature communications 


Article 


o) 


https://doi.org/10.1038/s41467-023-38523-4 


Hydrogen and dark oxygen drive microbial 
productivity in diverse groundwater 


ecosystems 


Received: 13 September 2022 S. Emil Ruff © 12° 


, Pauline Humez‘, Isabella Hrabe de Angelis ©", 


Muhe Diao ©', Michael Nightingale’, Sara Cho ©", Liam Connors ©', 


Accepted: 5 May 2023 


Published online: 13 June 2023 


® Check for updates 


Olukayode O. Kuloyo', Alan Seltzer ®°, Samuel Bowman®, Scott D. Wankel ©”, 
Cynthia N. McClain ©*®7, Bernhard Mayer ©' & Marc Strous ©' 


Around 50% of humankind relies on groundwater as a source of drinking 
water. Here we investigate the age, geochemistry, and microbiology of 138 
groundwater samples from 95 monitoring wells (<250 m depth) located in 14 
aquifers in Canada. The geochemistry and microbiology show consistent 
trends suggesting large-scale aerobic and anaerobic hydrogen, methane, 
nitrogen, and sulfur cycling carried out by diverse microbial communities. 
Older groundwaters, especially in aquifers with organic carbon-rich strata, 
contain on average more cells (up to 1.4 x 10’ mL") than younger ground- 
waters, challenging current estimates of subsurface cell abundances. We 
observe substantial concentrations of dissolved oxygen (0.52 + 0.12 mg L” 
[mean + SE]; n=57) in older groundwaters that seem to support aerobic 
metabolisms in subsurface ecosystems at an unprecedented scale. Metage- 


nomics, oxygen isotope analyses and mixing models indicate that dark oxygen 
is produced in situ via microbial dismutation. We show that ancient ground- 
waters sustain productive communities and highlight an overlooked oxygen 
source in present and past subsurface ecosystems of Earth. 


About 2% of Earth’s water resources occur as groundwater, of which 
half is saline and the other half fresh’. This subsurface freshwater, 
represents around 30% of the global freshwater resources, sixty times 
more than in all lakes, rivers, and the atmosphere combined, and 
exceeded only by the inaccessible and currently still frozen polar ice 
caps'. Aquifers and rock fractures may also hold up to 30% of the total 
microbial biomass on Earth’, contribute substantially to carbon 
fixation’, and contain high proportions of uncultured archaea, bac- 
teria, and viruses”> with a broad spectrum of lifestyles‘. Despite the 
global occurrence of groundwater and the magnitude and diversity of 
its resident biomass, our understanding of the composition and 


activity of microbial communities inhabiting these hidden aquatic 
ecosystems is still patchy, often developed from samples of a few 
select wells or a single aquifer. In particular, the geochemical and 
ecological processes that shape groundwater microbial communities 
over space and time are not well constrained’. 

Establishing robust links between microbial communities and 
geochemical characteristics of groundwaters requires large datasets 
having a comprehensive environmental inventory associated with 
each microbial community sample. To this end, the Groundwater 
Observation Well Network (GOWN) maintained by Alberta Environ- 
ment and Protected Areas (AEPA) in Canada, has compiled 
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geochemical data for over 250 groundwaters obtained from mon- 
itoring wells in different aquifers and geographic regions, representing 
a variety of geochemical regimes and groundwater ages. Each GOWN 
well has been sampled repeatedly for many years, including some for 
several decades®. Since 2006, this comprehensive monitoring pro- 
gram has systematically collected regular water level and chemical 
water quality information and isotopic compositions for aqueous and 
gaseous samples’. The province of Alberta is situated in the Western 
Canadian Sedimentary Basin, which hosts major oil, gas, coal, as well as 
sulfur, salt, limestone, and dolomite deposits. The shallow and deep 
subsurface has been extensively studied in the context of petroleum 
and coal exploration and development" (Fig. 1). 

In this work, we analyzed 138 groundwaters from 95 GOWN wells 
completed in bedrock and surficial aquifers across Alberta (Fig. 1, 
Table S1, Supplementary Data 1) with the aim of investigating the 
biogeochemistry and microbial ecology of a broad range of aquifer 
environments. We employed a multidisciplinary characterization of 
the groundwater geochemical composition including major and 
minor ion concentrations, dissolved gases, determination of 
groundwater ages, and the compositions and metabolic capabilities 
of its resident microbial communities. The guiding objectives were to 
determine (i) whether geochemistry is consistent with microbial 
community compositions and metabolisms, (ii) which electron 
donors and acceptors support these communities, and (iii) which 
microbial lineages represent key populations. We show that ancient 
groundwaters can harbor productive and diverse microbial com- 
munities. Metagenome analyses and geochemistry reveal that these 
communities make a living using hydrogen, methane, sulfur, and 
molecular oxygen. We show that molecular oxygen was likely 
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Fig. 1 | Sampling locations, geological formations, and groundwater ages. 
a Location of studied groundwater wells within the energy resources context of the 
province of Alberta. Colors indicate the groundwater age at each well (yellow: 
younger waters; red: intermediate age; blue: older waters that are sulfate-rich; 
purple: older waters with little sulfate). Circle size represents average microbial cell 
numbers in the groundwater samples, ranging from 10* (smallest full circle) to 10” 


produced in situ and plays a role in the biogeochemistry and ecology 
of these subsurface ecosystems. 


Results and discussion 

Geochemical evolution of groundwater in the Canadian Prairie 
The mean residence time of water in an aquifer is a very important 
factor influencing its geochemical composition (or facies) (Fig. 2a). 
Mean residence time, or groundwater age, is determined using 
radioisotopes such as tritium and “C and refers to the travel time 
between the point of infiltration and the point of sampling”. In this 
study, groundwater age inversely correlated with the ratio between 
calcium and sodium concentrations in groundwater (Fig. 2b). The 
youngest, tritium-containing groundwaters (indicating recharge after 
1952) were characterized by low total dissolved solids (<400 mg L”) 
and high Ca/Na ratios (median: 3.5). Dissolution of carbonates during 
infiltration led to high calcium, magnesium, and bicarbonate con- 
centrations in these waters. These young groundwaters (yellow sym- 
bols, Fig. 2) generally had low methane and sulfate concentrations 
(Fig. 3a, b) and were predominantly collected from wells completed in 
surficial Neogene-Quaternary deposits (Fig. 1b). 

In contrast, the oldest groundwaters (purple symbols, Fig. 2) 
contained no tritium and little “C, indicative of groundwaters more 
than several hundreds or thousands of years old. They had elevated 
total dissolved solids (>900 mg L”) and a low Ca/Na ratio (median: 
0.01). Old groundwaters were characterized by reducing conditions 
and contained high dissolved methane concentrations 
(12.8+2.4mgL7?  [mean+SE]; median: 0.72mgL7, range: 
0.001-74.2 mg L7; Fig. 3a, Supplementary Data 1). These waters had 
elevated sodium, bicarbonate, and chloride concentrations resulting 
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proportion of water types in the surficial, channel, and bedrock sediments, as well 
as in major geological formations of Alberta, showing that groundwater geo- 
chemistry evolved with the increasing age of the formations. NA not assessed, HSC 
Horseshoe canyon, Gp. Group. 
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Fig. 2 | Groundwater geochemistry and age dating. a Piper diagram of hydro- 
chemical facies of each groundwater sample (circles) visualizing calcium/sodium 
mass ratio used as a proxy for geochemical evolution, i.e., residence time/age. Na: 
sodium, K: potassium, Ca: calcium, Mg: magnesium, CI: chloride, S04”: sulfate, 
CO37: carbonate, HCO3: bicarbonate b Ca/Na ratio decreases with increasing 


from water-rock interactions, including ion exchange, and weathering 
of minerals. The older groundwater samples were obtained from wells 
completed in buried river valleys (channels) and Paleogene and Cre- 
taceous sedimentary bedrock formations that are often characterized 
by the presence of coal and/or shale”. 

Sulfate, predominantly derived from pyrite oxidation and to a 
lesser extent from anhydrite or gypsum dissolution”, was ubiquitous 
in a third group of groundwater samples (blue symbols, Fig. 2). These 
waters were characterized by long residence times, contained even 
higher dissolved solids (>1700 mg L”) and had intermediate Ca/Na 
ratios (median of 0.12). Sulfate was often the most abundant anion and 
electron acceptor in this group of groundwaters, resulting in sulfate- 
rich hydrochemical facies with low methane concentrations. These 
groundwater samples were collected from wells completed in surficial 
deposits, but also from bedrock aquifers completed in clastic, often 
marine sedimentary rocks of the Bearpaw formation (Fig. 1b). 

The mixing of groundwater of different ages and geochemical 
water types (Fig. 2a) often results in characteristic intermediate 
hydrochemical facies (red symbols, Fig. 2). These waters were mostly 
obtained from aquifers in surficial deposits and occasionally from 
bedrock aquifers (Fig. 1b). A detailed characterization of the geologic 
formations, groundwater geochemistry, and ages are provided in the 
Supplementary Results and Supplementary Data 1. 


Old groundwaters contain biomass-rich microbial communities 
Microbial cell density commonly decreases with increasing depth in 
marine” and terrestrial’ subsurface ecosystems. This decrease in cell 
numbers and biomass is often attributed to energy limitation in sub- 
surface ecosystems". Within aquifers cell numbers may not necessarily 
decrease with increasing depth”, but due to other factors including 
increasing water age’. To our surprise, we found significantly more 
cells in aquifers with old groundwater than in those with younger water 
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residence time (“*Cpic uncorrected age, BP: before present). “C: “carbon, DIC: 
dissolved inorganic carbon. Samples that accumulate in the plot around 40k yrs BP 
(detection limit of method) may be much older. 7H-positive samples (circles with 
crosses) comprise meteoric water past 1952 (H-bomb) and corroborate the age 
trend. c Schematic timeline and summary of water aging. 


(Fig. 4, Supplementary Data 2). This suggests that the older aquifers 
and geochemically evolved groundwaters are in fact productive eco- 
systems that provide energy for the growth of microorganisms. 
Average cell numbers in the geochemically evolved groundwater 
samples reached 10’ cells mL” (Fig. 4). Overall, microbial cell numbers 
across the entire studied region, did not decrease with depth (Sup- 
plementary Fig. 2a, b). Cell counts were on average slightly higher in 
aquifers in geologic strata that contained shales or coal (Supplemen- 
tary Fig. 2c), suggesting that elevated contents of organic carbon may 
provide additional energy sources to microbes. It should be noted that 
the cell numbers presented here are conservative estimates likely 
capturing only the free-living cells from water samples which cannot 
account for the substantial number of cells living in biofilms”°. 

The morphology and size of cells varied greatly, and included 
cocci-, rod-, vibrio-, and spiral-shaped cells, as well as filaments and 
small aggregates (Supplementary Fig. 3). Cell morphologies, sharp cell 
boundaries and the bright signal from nucleic acid staining suggested 
a largely active community”. Average cell size was similar across 
groundwater samples and hence high cell numbers are a proxy for high 
biomass and productivity (Supplementary Fig. 3). To the best of our 
knowledge this is the first time that consistently high cell abundances 
were documented in old aquifers across a large geographic area 
(>210,000 km”). 

The increase in cell numbers with groundwater age was accom- 
panied by a substantial decrease in archaeal and bacterial diversity 
(Fig. 5a, b; Supplementary Fig. 4) and a shift in microbial community 
structure from young to old groundwaters (Fig. 5c, d; Supplementary 
Fig. 5). Of the environmental parameters that we tested, the ones that 
most strongly explained the variance in microbial community struc- 
ture were methane concentrations and geochemical proxies for 
groundwater age, including the concentrations of sodium, calcium, 
and magnesium and well depth (Fig. Se, f). Decreases in diversity and 
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Fig. 3 | Groundwater isotope geochemistry. Microbes mediating methane (CH4) 
and sulfate (S04%) cycling impact the carbon and sulfur pools in the groundwater 
systems. Each circle represents a sample. Circle color depicts water age. Trends 

concerning reduction and oxidation of the compounds are indicated by arrows or 
areas. a Carbon isotopic signature of methane versus methane concentration. PDB 
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shifts in community structure are common features of microbial 
bloom situations in productive ecosystems, which often feature ele- 
vated abundances of a few key species””. 


Hydrogen as basal energy source in productive aquifers 
Our results indicated that a major energy source driving high pro- 
ductivity in mature groundwater is hydrogen. Considerable amounts 
of dissolved hydrogen were detected at selected sites (0.002, 0.007, 
and 0.1% in samples GW121, GW223, and GW951), and in many aquifers 
hydrogenotrophic methanogens affiliating with Methanobacterium, 
Methanoregula, and Methanospirillum showed high relative abun- 
dances of 16S rRNA gene amplicon sequence variants (ASVs; Fig. 6a, 
Supplementary Figs. 6, 7, Supplementary Data 3). Hydrogen serves as 
sole electron donor for these lineages and the isotope composition of 
CH, and CO, also supports widespread hydrogenotrophic methano- 
genesis (Fig. 3c). Methane with concentrations ranging from 0.006 to 
more than 74 mg L™ was present in all analyzed groundwater samples 
(Fig. 3a). Most samples, however, had concentrations <1 mg L” (73%; 
n=106, Supplementary Fig. 6). Based on combined carbon and 
hydrogen isotope composition, the low methane concentrations likely 
reflect high relative consumption by microbial methane oxidi- 
zers (Fig. 3a). 

Among the detected methane oxidizers were obligate anaerobic 
methane-oxidizing archaea (ANME) of the genus Methanoperedens 


(ANME-2d, Fig. 6a). These methanotrophs were often found in aquifers 
that contained either dissolved iron or manganese or both. ANME-2d 
were found in more than half of the archaeal ASV datasets (34 of 64), 
predominantly in methane-rich and sulfate-poor older groundwaters, 
in line with their previously observed occurrence in deep granitic 
water environments”. Their ability to oxidize methane using iron and 
manganese oxides” could explain the high relative ASV abundance of 
over 40% in the sulfate-depleted mature groundwaters of at least four 
samples (GWI111, GW218, GW311, and GW456). Methanoperedens was 
the fifth most abundant archaeal genus in the studied aquifers and has 
the metabolic potential to connect the carbon, nitrogen, and metal 
cycles via the anaerobic oxidation of methane”*””. Anaerobic oxidation 
of methane coupled to bacterial sulfate reduction might occur in 
certain aquifers as well, because ANME-2ab made up almost half of the 
archaeal community at GW160, and ANME-2c were detected at GW220 
(Supplementary Data 3). Bacterial Methanomirabiliales, an order 
known to comprise microbes that are capable of oxidizing methane in 
anoxic environments using nitrite’, were abundant in many aquifer 
datasets, especially in young and intermediate waters (Supplemen- 
tary Fig. 8). 

Hydrogen also could have sustained the growth of organisms 
affiliating with the genera Desulforudis and Desulfomicrobium and 
many other sulfate and sulfur reducers, and sulfur disproportionators, 
which were diverse and widespread in the studied groundwater 
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Fig. 4 | Cell abundances in groundwaters. Cell abundance was determined using 
fluorescence microscopy and is given on a logarithmic scale. Boxplots show 1* and 
3" quartile, maximum and minimum values (whiskers), median (line), mean (tri- 
angle), and outlier (dots). a Boxplots represent groundwater samples from indivi- 
dual wells and summarize cell counts from n independent fields of view. n = 40: 
Samples 3-10, 12, 15-22, 24-29, 34-47, 49-51, 54, 56, 63, 67-68, 70-71. n = 39: 
Samples 1, 14, 33. n = 38: Sample 2. n = 20: Samples 13, 31-32, 55, 57, 61, 72. n=12: 
Sample 73. For more details on individual cell counts see Supplementary Data 2. X 


axis: Numbers correspond to well IDs (due to space constraints well IDs are given in 
Supplementary Data 1). Wells which were sampled at two different time points are 
marked with a star (*), technical replicates are marked with a plus (+). b Boxplots 
represent cell numbers averaged across groundwater age summarizing the average 
cell numbers in a (excluding technical replicates). n number of samples. Sig- 
nificance was tested using a Wilcoxon rank sum test. Significance levels are: 

*p < 0.05; **p < 0.01; **p < 0.001; (uncorrected). 


ecosystems (Fig. 6b, Supplementary Fig. 9, Supplementary Data 4). 
Candidatus Desulforudis audaxviator were the most abundant sulfate 
reducers across all ASV datasets, yet predominantly occurred in old 
groundwaters. They often made up >50% of the clades known to 
comprise sulfur-cycling microbes. Ca. Desulforudis are hydrogen-oxi- 
dizing, sulfate-reducing Clostridia reported to thrive in deep terrestrial 
aquifer ecosystems”. In samples GW3026 and GW217 we found high 
sequence abundances of microbes affiliating with obligately syn- 
trophic Smithella sp. and Syntrophus sp. that live together with 
organisms scavenging hydrogen*””. 

Shotgun metagenomics supported the inferences derived from 
amplicon sequencing and revealed genes encoding Ni-Fe hydro- 
genases, hydrogenotrophic methanogenesis and sulfate reduction in 
the groundwater communities (Fig. 7). A high-quality metagenome- 
assembled genome (MAG) of a population within the Methanobacter- 
iaceae (UBA349, MAG-32) contained a complete methyl-coenzyme M 
reductase (mcrABCG) and all other genes needed for hydrogenotrophic 
methanogenesis (cdhA-E, ftr, mch, fae), yet lacked genes diagnostic for 
methylotrophic, or aceticlastic methanogenesis (Supplementary 
Data 5). A complete sulfate reduction pathway was found in MAG-18 an 
organism in the genus-level clade UBA2262 within the Desulfovi- 
brionaceae (Fig. 7). The hydrogen could be derived from microbial 
fermentation, e.g., via Fe-Fe hydrogenases that were present in three 
MAGs affiliating with Bacteroidia (Lutibacter, UBA6024, JAAYJTO1; 
Supplementary Data 5). These MAGs also encoded the TonB/SusC 
transport system, which many Bacteroidia use to consume oligo- 
saccharides. Hydrogen can also derive from abiotic sources including 
shales, coal beds, and other organic-rich strata”, water-rock reactions”, 
pumping equipment or steel well casings”, or even water radiolysis*”*. 


Aquifers harbor microbes affiliating with aerobic and faculta- 
tively anaerobic lineages 

We found high relative ASV abundances of microbes affiliating 
with obligate and facultatively aerobic hydrogen-, methane-, and 


sulfur-oxidizing lineages in most of the old groundwaters obtained 
from confined aquifers (Figs. 6, 8, Supplementary Figs. 7-10, Sup- 
plementary Data 4). The presence of potentially aerobic organisms 
was not the result of sample handling, contamination, or of microbial 
growth after sampling, and alpha and beta diversity, composition, 
and cell counts did not change with sample storage time (Supple- 
mentary Fig. 11, Methods). Thus we are confident that these microbes 
were present in the aquifers at the time of sampling. Analyses of 
metagenomes from five wells (GW114, GW144, GW218, GW265, and 
GW972, all bearing old groundwaters) corroborated the community 
composition detected via metabarcoding. Mapped metagenomic 
16S rRNA gene short reads (Supplementary Fig. 12), as well as 
reconstructed full-length 16S rRNA gene sequences (Supplementary 
Data 6) showed the same key genera and often the same species as 
the ASV datasets. The reliability of the ASV-based community com- 
position was further supported by very high sequence similarities 
between ASVs and reconstructed full-length 16S rRNA genes (Sup- 
plementary Fig. 13, Supplementary Data 7). Both sequencing 
approaches show the same key players, including Methylobacter 
tundripaludum (MAG-01 - 04), Methylotenera sp. (MAG-09 - 13), 
Hydrogenophaga sp. (MAG-14 - 17), Sulfuricurvuum sp. (MAG-19) and 
Thiobacillus sp. (MAG-20; Fig. 7). 

We detected Hydrogenophaga at high relative abundances in most 
aquifers (Figs. 6b, 7, 8a). Hydrogenophaga MAGs contained genes 
encoding hydrogenases, oxidases, and RuBisCO, but also a sulfur 
oxidation pathway (sqr, fccAB, soeABC, soxABX), a complete deni- 
trification pathway (in 2 of 4 MAGs) and an aerobic carbon monoxide 
dehydrogenase (coxSML), suggesting a facultatively anaerobic 
lithoautotrophic lifestyle. Microorganisms affiliating with aerobic or 
denitrifying methano- and methylotrophic bacteria” were also abun- 
dant in all aquifers, especially in older groundwaters (Fig. 6b, Supple- 
mentary Fig. 8). Methylotrophic Methylotenera sp. were the most 
abundant populations in the ASV datasets, followed by methane- 
oxidizing members of the genus Methylobacter. Both genera were also 
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Fig. 6 | Microbial community composition. Relative sequence abundance of the 
top 20 most abundant a archaeal and b, c bacterial lineages at genus level based on 
16S rRNA gene amplicon sequence variants. The remaining lineages with lower 

abundances are summarized as Other. Genera that were unclassified in the SILVA 
reference database are abbreviated as unc and the next higher phylogenetic level is 


@e 
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shown. The two archaeal clades denoted by * belong to the order Woesearchaeales. 
One methanogen, as well as almost half (9 of 20) of the top bacterial genera are also 
represented by metagenome-assembled genomes (Fig. 7). The legend in b applies 
to all panels. 


abundant in the metagenomes of the five wells and were represented 
by nine MAGs. Methylobacter MAGs had complete pathways for the 
oxidation of methane to carbon dioxide including the critical activa- 
tion of methane with molecular oxygen via the particulate methane 
monooxygenase (pmoABC). It was recently suggested that Methylo- 
bacter can couple methane oxidation to denitrification in hypoxic 
conditions”. Indeed, we found these organisms to be capable of partial 
denitrification, from nitrite to nitrous oxide. Methylotenera MAGs 
contained genes to oxidize methanol to carbon dioxide and reduce 
nitrate to nitrite (Fig. 7). Denitrification by Methylotenera was reported 


previously” and hence Methylotenera could provide nitrite to Methy- 
lobacter, while Methylobacter provides methanol to Methylotenera. 
Complementary interactions between Methylobacter and Methylote- 
nera were observed in a methane-rich shallow aquifer*’, a lake“ as well 
as in methanotrophic enrichment cultures*’. The metabolic cap- 
abilities of the consortium however suggest that molecular oxygen 
needs to be available to catalyze the initial activation of methane, even 
if subsequent steps can be coupled to denitrification. Methyli- 
corpusculum sp. (MAG-05), in contrast seemed to be entirely lacking 
the capability to use other electron acceptors than oxygen (Fig. 7). 
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Fig. 7 |Metagenome-assembled genomes and key metabolic pathways. Selected 
metagenome-assembled genomes (MAGs) of abundant microbial populations in 
groundwaters show the genetic capabilities of methane, sulfur, and nitrogen 
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Aerobic methane oxidizers affiliating with the genus Crenothrix*® were 
the fourth most abundant methanotroph and occurred in almost half 
of the ASV datasets (45 of 109, Supplementary Fig. 8, Supplementary 
Data 6). Crenothrix sp. are filamentous microbes that can occur in 
groundwater“ and may represent the filamentous cells observed with 
the microscope in many samples (Supplementary Fig. 3). We also 
retrieved three high-quality MAGs of the uncultured genus JABFRCO1 
within the Methylomonadaceae (MAG-06 - 08). These MAGs encoded a 
soluble and a particulate methane monooxygenase and complete 
methane oxidation pathway revealing that the organism is an aerobic 
methane oxidizer requiring oxygen for the initial oxidation of methane 
to methanol. Like the closely related Methylobacter the organisms 


affiliating with JABFRCO1 may use nitrite as an alternative electron 
acceptor for the further oxidation of methanol to carbon dioxide in 
hypoxic conditions (Fig. 7, Supplementary Data 5). Of the twenty 
most abundant bacterial genera known to oxidize one-carbon com- 
pounds (Supplementary Fig. 8), twelve were associated with oxida- 
tion of methanol, methylamines, and other methyl compounds“, 
while eight were associated with methane oxidation*®. The wide- 
spread presence of microorganisms that are capable of aerobic 
methanotrophy corroborates previous findings of aerobic metha- 
notrophic activity in coalbed formations*’**, with our study 
expanding this finding to aquifers in more diverse geological settings 
and a much larger geographical scale. 


Nature Communications | (2023)14:3194 


Article 


https://doi.org/10.1038/s41467-023-38523-4 


100 


50 


Relative sequence abundance 
of Hydrogenophaga (in %) 


100 | 


a 
fo) 


Relative sequence abundance 
of Methylomonadaceae (in %) 


Relative sequence abundance | 
of aerobic/denitrifying sulfur oxidizers (in %) 
a > 
7 S 


fo) 


I Younger 
E a I Internediate 
i im ° * | I Older w/ SO, 


— 


0.1 1 10 
Dissolved oxygen (mg L”) 


E Older 
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trifying clades tend to peak at hypoxic conditions (-0.5 mg L”) potentially 
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Sulfur-oxidizing bacteria including Sulfuricurvum sp., Thiobacillus 
sp., Thiomicrorhabdus, and Sulfurimonas were also widespread and 
abundant in the studied aquifers (Fig. 6c, Supplementary Fig. 10). 
Thiobacillus (MAG-20) was able to oxidize sulfide and thiosulfate to 
sulfate using oxygen, but also had a complete nitrate reduction path- 
way (Fig. 7). A similar lifestyle was likely for PFJXO1 within the Thio- 
bacillaceae (MAG-21), Sulfurimicrobium (MAG-22), Rhodoferax (MAG- 
23, -24), UBA2250 (MAG-26, —27) and ETT8 (MAG-29) affiliating with 
Rhodocyclaceae and Rhodobacteraceae. Some members of these 
clades were shown to be active facultatively anaerobic denitrifying 
sulfur oxidizers in the terrestrial subsurface**°°. The versatile sulfur 
metabolisms of these clades could explain previous observations of 
pyrite oxidation, as well as the presence of sulfur reducers like 
Desulfuromonas which can reduce elemental sulfur” and replenish the 
sulfide pool. 

We found obligate aerobic ammonium oxidizing Nitrosomonas, its 
MAG encoded a partial pmoA/amoA and a hao to oxidize ammonium to 
nitric oxide (Supplementary Data 5). Like the activation of methane, 
the ammonium monooxygenase requires molecular oxygen. The bio- 
mass produced by the diverse, coexisting autotrophs in turn may 
support multitudes of archaeal and bacterial heterotrophs, including 
DPANN archaea and Patescibacteria. Numerous MAGs represented 
Flavobacteriales, Cytophagales, Chitinophagales, and Bacteroidales, 
heterotrophs with an appetite for oligosaccharides and the necessary 
transport systems. In the archaeal community aerobic Woesearch- 
aeales were particularly diverse and abundant in all aquifers except 
those with old groundwaters, and in the bacterial community poten- 
tially heterotrophic Comamonadaceae and Pseudomonadaceae were 
also very abundant (Supplementary Fig. 14b, c). Heterotrophic clades 
remineralize organic carbon compounds to carbon dioxide, which 
could then be recycled by hydrogenotrophic methanogens and other 
autotrophs to produce biomass. 


Dark oxygen is an electron acceptor in aquifers 

As water infiltrates through recharge zones and flows through the 
aquifers, the dissolved oxygen (DO; assumed to be at saturation 
equilibrium initially) is consumed by microbial respiration, decom- 
position of organic matter, or by reacting with reduced minerals”. 
These oxygen-consuming reactions often reduce the DO content of 
groundwater to below detection limits, particularly in groundwater 
with long residence times that has been out of contact with the 
atmosphere for many years, centuries, or even millennia”. To our 
surprise we detected low concentrations of DO in most groundwater 
samples including those of deeper aquifers containing geochemically 
mature groundwater (Fig. 9a). Using the known stoichiometry of 
relevant microbial metabolisms, we estimated that the consumed 
oxygen (DO saturation at infiltration minus measured DO concentra- 
tion in the obtained samples) could sustain more microbial cells than 
we observed in 85% of aquifers (59 of 70; Supplementary Data 8). In the 
remaining 11 aquifers, the DO content was below the detection limit. 
The occurrence of more than 0.3 mg L” DO in many old groundwaters 
of deeper and confined aquifers was remarkable and suggests that 
oxygen is present in ecosystems that are often assumed to be 
anoxic’. 

To investigate whether oxygen may have been introduced during 
sampling, we carried out oxygen isotope analyses of molecular oxygen 
(O2). Indeed, some oxygen isotope data for DO were consistent with 
groundwater that is in equilibrium with atmospheric oxygen 
(600, =~+23.9%o + 0.1%o, Fig. 9b), indicating air contamination during 
sampling, or the presence of atmospheric-derived oxygen in the 
aquifers. However, at certain sites (GW218, GW265) we observed 
markedly lower "Oo, values (as low as +21%o), while also finding 
elevated O,:Ar ratios (Fig. 9b). The lower 5°00, and higher O,:Ar ratios 
interestingly fell along a trend, representing the simulated addition of 
DO with a "Oo, value much lower than that of air-equilibrated water 
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Fig. 9 | Concentration and isotopic signature of oxygen in groundwater sam- 
ples. a Depth profile and dissolved oxygen concentration (mg L”) in the ground- 
water samples. b 180-0, isotopic signature over oxygen:argon ratio (O/Ar). The 
composition of lab air-equilibrated water and lab air are represented by a gray and 
black circle, respectively. Error bars are 1 standard deviation (based on repeat 
analysis of lab air, n = 13) and are smaller than the symbols. The blue and black solid 
lines show different isotope effects that would be caused by fractionation (e) due to 
respiration (resp.). During consumption of O2 by respiration (decreasing O2:Ar 


(dashed trend line in Fig. 9b). The variation among the triplicate sam- 
ples may be explained by the variable presence of atmospheric-derived 
oxygen and/or by varying microbial respiration, which will decrease the 
O,:Ar and increase 5*0o, values. In fact, mixing and biological pro- 
cesses are likely both at work. Most remarkably, however, there is no 
plausible scenario to our knowledge whereby air contamination or 
microbial respiration could lead to the observed increasing trend in 
O,:Ar ratio coupled with decreasing 5°Oo, values. That is, air con- 
tamination would bring O2:Ar ratios and 5°0o, back toward air- 
equilibrated values of +23.9%o° while microbial consumption of oxygen 
would lower O;:Ar while increasing 5'°O,,°°. The most parsimonious 
explanation for the observed trend would be the in-situ production of 
O, with a very low 5*0o, (e.g., -20%o). Given the low &'80 (of H20) of 
groundwaters in Alberta (—18.6%o + -2.0%o; mean + SD; n=144, Sup- 
plementary Data 1), which would be the source of oxygen for any in-situ 
production of O, even a small amount of production would have a 
large impact on lowering the overall 5*0o, in groundwater, in the 
direction of our observations. In the absence of light, oxygen can be 
produced by water radiolysis”, or microbially by chlorite dismutation’®, 
nitric oxide dismutation’®”’, and H202 dismutation®’. 


Microbial production of dark oxygen via dismutation 

Microbes can produce O; via chlorite dismutation, a process carried 
out among others by microbes of the genera Dechloromonas*, 
Dechlorobacter, Dechlorosoma, Azospira, Azospirillum’, Nitrospina and 
Nitrobacter, all of which were present in the GOWN aquifers based on 
the 16S rRNA gene analyses. ASVs affiliating with the genus Dechlor- 
omonas were among the most abundant in the study accounting for up 
to 30% of all bacterial sequences in some samples. We found 13 chlorite 
dismutase genes mainly from Dechloromonas, Nitrospina, Nakamur- 
ella, and Magnetospirillum species in old, anoxic/hypoxic (DO: 
0-0.55 mg L”) groundwaters of GW144, GW218, GW265, and GW972 
(Supplementary Data 9). The chlorite dismutase genes were detected 
at considerable sequencing depths ranging from 4 to 212x, indicating a 
moderate to high abundance of the organisms having this gene. Par- 
ticularly two chlorite dismutases assigned to Magnetospirillum had 
very high sequencing depth (207x and 212x). Groundwater can also 
contain a high diversity of nitric oxide dismutation genes” and 
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ratio), there is a preferential accumulation of O in the remaining O- pool—leading 
to higher 5°05, values. As the degree of this isotope fractionation can vary the two 
solid lines represents two hypothetical but realistic scenarios for how the O pool 
might change with microbial respiration. The black dashed line represents a mixing 
line between air-equilibrated water mixed with a hypothetical ‘new source of Oy 
that has a very low "Oo, (-20%o vs V-SMOW). Such isotopically light oxygen is 
consistent with the biological formation of O2. V-SMOW Vienna Standard Mean 
Ocean Water. 


microbes affiliating with Methylomirabilaceae that can potentially 
dismutate nitric oxide were present in many samples. We also found 
nitric oxide dismutase genes (nod) in high-quality MAGs of Sedimini- 
bacterium (MAG-30) and Algoriphagus (MAG-31; Fig. 7; Supplementary 
Data 9). Both lineages were recently shown to have a nod®. The nod 
affiliating with Sediminibacterium was sequenced at a depth of 700x 
corroborating the very high abundance of Sediminibacterium in the 
community (12%). Nitric oxide reductase is related to nod but it can be 
identified and distinguished by several diagnostic amino acid residues 
in the active center of the enzyme (Supplementary Fig. 15). It was 
recently shown that nitric oxide dismutating Nitrosopumilus sp. can 
leak intracellularly produced oxygen into the medium possibly sup- 
porting other aerobic organisms”. It was also reported that even 
Pseudomonas aeruginosa seem to be able to dismutate nitric oxide and 
release peaks of up to 20 uM oxygen into their surroundings™. Linea- 
ges affiliating with both, Nitrosopumilus and Pseudomonas were wide- 
spread and abundant in the studied groundwaters. Abundant 
membrane-bound dismutases in the studied aquifers and their 
release of oxygen into the surroundings could explain the observed 
presence of O-depleted oxygen in old groundwaters. 


Microbial lineages, metabolic capabilities, and niches in 
groundwater ecosystems 

Our findings suggest that hydrogen is a basal energy source fueling a 
rich mosaic of microbial metabolisms. The extensive microbial pro- 
ductivity in the studied groundwater ecosystems (Fig. 10) underlines 
the importance of hydrogen in the terrestrial subsurface®. Hydro- 
genotrophic methanogens produce methane, which in turn explains 
the high abundance of methanotrophs and widespread genes involved 
in methanotrophy. Substantial aerobic and anaerobic methane oxi- 
dation can reduce greenhouse gas emissions from aquifers containing 
coal beds and shales with relevance for carbon budgets and climate 
change. Hydrogen also appears to drive extensive microbial sulfur 
cycling. Sulfate reducers produce sulfide which is then oxidized by 
thiotrophic microbes using either oxygen or nitrate as an electron 
acceptor in their metabolic pathways. The microbial community 
composition, the detected metabolic capabilities, and aqueous and 
isotope geochemistry are all aligned and suggest that groundwater 
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microbial communities utilize the entire breadth of known biochemi- 
cally available redox potential, from hydrogen to oxygen. Unlike in 
many other aquifers in which cell numbers decrease with increasing 
groundwater age, cell numbers increased with age in groundwater 
obtained from aquifers in Alberta, corroborating that these aquifer 
communities are fueled by autochthonous energy sources derived 
from within the subsurface. Considering the size of the investigated 
area we conclude that global subsurface biomass may be under- 
estimated, particularly in aquifers containing organic carbon-rich 
strata, such as coal beds and shales. The geochemical and micro- 
biological data suggest that oxygen is an important electron acceptor 
for microbial metabolisms in the studied aquifers. The presence of 
aerobic microbes in deep aquifers and bedrock has been previously 
reported®®, yet the mechanisms of deep oxygen migration or in situ 
production remain unclear. While co-migration of some oxygen with 
the aging groundwater could theoretically explain the observed 
microbial productivity, the presented oxygen isotope analyses, and 
the abundance of microbial dismutase genes capable of producing 
oxygen suggest that at least some portion of the oxygen is generated 
in-situ. Most of the hypoxic waters surveyed here are several thousand 


to >10,000 years old based on tritium and “C data, supporting that 
even deep subsurface ecosystems provide niches for aerobic 
microorganisms”. The production of dark oxygen that we postulate in 
this work could provide a mechanism for previously reported iso- 
topically light groundwater oxygen anomalies that have lacked a clear 
explanation thus far?**°*”, Microbial dark oxygen production in 
subsurface ecosystems may thus be relevant for the functioning and 
evolution of the geobiosphere, as it provides a source of oxy- 
gen independent of light, on Earth as well as potentially on other 
celestial bodies. 


Methods 

Groundwater wells 

Approximately 40-90 groundwater wells are sampled annually for 
water quality by Alberta Environment and Protected Areas? and ana- 
lyzed for a variety of geochemical, isotopic, and microbial parameters. 
For this study, we analyzed 138 groundwater samples from 95 mon- 
itoring wells, sampled from January 2016 to July 2020. 38 wells are 
screened within Neogene-Quaternary surficial deposits while 57 wells 
reach sedimentary bedrock, usually of Paleogene or Cretaceous age 
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(Supplementary Data 1). Well depths in this study ranged from five to 
231m (mean=60m). The sedimentary bedrocks in which ground- 
water monitoring wells are completed include the Paskapoo and 
equivalent (n=25), Horseshoe Canyon and equivalent (n =13), Belly 
River Group (n= 9), Bearpaw (n= 4), Milk River (n = 3) and Loon River 
(n=2) formations. One well was not assigned with a geological for- 
mation. For further details concerning the location and geology refer 
to the SI. 


Sampling 

All wells were sampled at least once, while select wells were sampled 
during different seasons (up to three), or repeatedly to obtain sam- 
pling replicates (up to three). Some wells were installed as clusters, 
with multiple wells (up to three) at the same location completed at 
different depths within the same formation or across multiple forma- 
tions. All wells were purged to improve sample quality” and samples 
were collected after field parameters (dissolved oxygen - DO, 
oxidation-reduction potential - ORP, temperature, electrical con- 
ductivity - EC) had stabilized, indicating representative aquifer water. 
All wells were sampled following the same procedure and the samples 
were stored at 4 °C until further processing. Dependent on the location 
of the wells, the samples were in transit at 4 °C between 1 and 7 days. 
Tritium CH), and carbon-14 of dissolved inorganic carbon (“Cpyc) 
samples were collected with a submersible pump (SP). Samples for 
other water quality parameters, isotopes, dissolved gas, and micro- 
biological samples were collected with a bladder pump (BP). The SP is 
placed 2 m below the expected drawdown level for purging and sam- 
ple collection. The submersible pump has polyethylene tubing that 
does not get cleaned, and Nalgene tubing that has been cleaned with a 
general laboratory detergent (Sparclean) followed by a rinse with 
ultrapure (18 MQ) water. The BP is either placed in the vicinity of the SP 
or just above the screened interval for sample collection and has acid- 
washed, methanol, and ultrapure water-rinsed Teflon tubing. For 
microbiological analyses we collected 1000 mL” groundwater in a 
sterile Nalgene bottle. The sample was sent to the laboratory at 4 °C 
and used for microbial cell cryo-preservation, microbial cell counts, 
and nucleic acid extractions. 


Dissolved oxygen measurements 

Dissolved oxygen measurements were conducted with multi- 
parameter continuous water quality monitoring sondes, which have 
lifetime-based optical luminescence sensors with a resolution of 
0.01 mg L” and an accuracy of +0.1 mg L”. The specific sondes and 
sensors used during the 2016-2017 sampling campaigns were a YSI 
6600 XLM (Dissolved Oxygen - ROX; YSI, Yellow Springs, USA), EXO1 
or EXO2 (Optical Dissolved Oxygen Smart Sensor; YSI, Yellow Springs, 
USA), In-Situ AquaTROLL 600 (Rugged Dissolved Oxygen Sensor with 
RDO-X cap; In-situ, Fort Collins, USA), or an Ott Hydrolab DS5X (Hach 
LDO; Hach, Düsseldorf, Germany). The sample tubing was connected 
to a -400-900 mL flow through cell and field parameters were logged 
at 1-2 min intervals such that a full volume of water was displaced in 
between readings, using the default or accelerated averaging modes 
(5-40 second filtering), until parameter stabilization (+0.2-0.3 mg L”) 
was reached for 15 consecutive readings. 


Dissolved oxygen sensor calibration 

Before each field trip, calibrations of the dissolved oxygen sensors 
were performed in the lab. Calibrations for water-saturated air were 
conducted according to the user manuals of the probes YSI EXO or 
6-Series (YSI, Yellow Springs, USA) and the Aqua TROLL 600 (In-situ, 
Fort Collins, USA) or air-saturated water according to the user manual 
of the Hydrolab DSX6, DS6 and MS6 (Hach, Dusseldorf, Germany). For 
water-saturated air, the calibration cup was filled with 1/8 to 1⁄2 of an 
inch of tap water. A clean and dry dissolved oxygen sensor was 
inspected for damages and cracks, then covered with a clean 


calibration guard. The sensor and guard were inserted into the cali- 
bration cup and left loosely affixed to the probe to allow for venting. 
The setup was left to rest for 2-15 min out of direct sunlight to allow 
the temperature and pressure to equilibrate with water-saturated air. 
For air-saturated water, 1 L of room temperature water that had been at 
equilibrium with atmospheric pressure for at least 12 h, was vigorously 
shaken for one minute, poured into the calibration cup with the Hach 
LDO sensor cap (Hach, Düsseldorf, Germany) and the temperature 
sensor fully submersed. The calibration cup cap was lightly placed over 
the calibration cup to allow for air exchange and kept out of direct 
sunlight for 3-5 min. Barometric pressure was measured with a bar- 
ometer in an EXO Handheld (YSI, Yellow Springs, USA) or Vaisala 
PTB220 Digital Barometer (Vaisala Oyj, Vantaa, Finland) was recorded 
on the calibration record, and entered to the calibration software in 
mmHg. A one-point calibration was conducted for water-saturated air 
or air-saturated water. Readings were observed until they stabilized for 
40 seconds, the current value was recorded in the calibration record, 
and the calibration point was accepted. Percent saturation was con- 
verted to mg L” using measured temperature (thermistor) and salinity 
(conductivity sensor) and formulas from Standard Methods for the 
Examination of Water and Wastewater”. The health of the sensor cap 
was tracked by monitoring optical dissolved oxygen gain within 0.7 to 
1.4 and replaced every ~2 years. 


Chemical analyses 

Groundwater samples for major cations and anions were collected in 
polyethylene bottles and analyzed at ALS Global by Inductively Cou- 
pled Plasma-Mass Spectrometry (ICP-MS) and lon Chromatography, 
respectively. Total and phenolphthalein alkalinity was determined by 
titration, and bicarbonate and carbonate concentrations were subse- 
quently calculated. Dissolved organic carbon was measured by high- 
temperature combustion on a total organic carbon analyzer with 
infrared detection. Dissolved trace element samples were filtered in 
the field using 0.45 um in-line capsule filters, preserved to pH <2 with 
nitric acid, and analyzed at InnoTech Alberta by ICP-MS. 


Gas analyses 

Dissolved gas samples were collected with a bladder pump. Evacuated 
crimp-top glass serum bottles with butyl septa and preserved with 
mercuric chloride were filled by piercing the septa with a needle. Gas 
composition was determined in the Applied Geochemistry Laboratory 
at the University of Calgary by gas chromatography. The gas dryness 
parameter is defined as the ratio between the concentrations of 
methane and those of higher n-alkanes. The isotopic composition of 
methane and carbon dioxide was analyzed in the Isotope Science 
Laboratory at the University of Calgary on a ThermoFischer MAT 253 
isotope ratio mass spectrometer (IRMS) coupled with a Trace GC Ultra 
and GC Isolink (ThermoFisher) and reported relative to V-PDB for 6°C 
and V-SMOW for 67H. The precision for carbon isotope analyses was 
higher than +0.5%o for hydrocarbons, +0.2%o for carbon dioxide, and 
+2%o for &’H of methane. 


Groundwater dating 

All samples for “Cpic analyses were shipped to the A.E. Lalonde 
Laboratory in Ottawa for analysis by accelerator mass spectrometry 
(AMS). Sample pretreatment techniques can be found in ref. 74. 
Radiocarbon analyses are performed on a 3MV tandem accelerator 
mass spectrometer built by High Voltage Engineering (HVE). We used 
the fractionation-corrected fraction modern (the F“C value) for post- 
bomb “C data according to the amended conventions”. Tritium ÊH) is 
a radioactive isotope of hydrogen with a half-life of 12.4 years. Tritium 
concentrations are measured in tritium units (TU) where 1 TU is 
defined as the presence of one tritium in 10% atoms of hydrogen. For 
samples from the 2016-2017 sampling campaign, tritium in the 
obtained groundwater samples was analyzed at the A.E. Lalonde 
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Laboratory in Ottawa by electrolytic enrichment and the standard 
method of liquid scintillation counting with a precision of 0.8 TU. 


Carbon, oxygen, nitrogen, and sulfur isotope analyses 

Most stable isotope analyses were performed in the Isotope Science 
Laboratory at the University of Calgary (https://www.ucalgary.ca/labs/ 
isotope-science-lab/techniques). Dissolved inorganic carbon (DIC) 
isotope samples were filtered in the field using an in-line 0.2 um cap- 
sule filter and collected in a 125 mL clear glass flint bottle with a cone 
cap. CO, was produced from DIC by acid attack and &°C was analyzed 
as described in the gas analyses section. For the isotopic composition 
of sulfate, samples were filtered in the field with an in-line 0.2 ym 
capsule filter and collected in 1L wide mouth HDPE plastic bottles. 
Dissolved sulfate was converted into barium sulfate and analyzed 
using a ThermoQuest Finnigan Delta Plus XL IRMS coupled with a 
Fisons NA 1500 Elemental analyzer for 6°*Sso4 and a HEKAtech HT 
Oxygen Analyzer with a Zero Blank autosampler for 5'°Oso4. 8**Sso4 is 
reported relative to V-CDT and 8Oso4 is reported relative to V-SMOW. 
Precisions for both 88Oco4 and &’*Sso4 is +0.5%o. The isotopic com- 
position of nitrate was determined on N20 generated by the denitrifier 
technique using a Thermo Scientific Delta V Plus IRMS coupled with a 
Finnigan MAT PreCon. Precisions for 8°Nyo3 and "Oyos are +0.3%o 
and +0.7%o, respectively. 


Oxygen isotope analyses on dissolved oxygen 

Samples for the analysis of oxygen isotope ratios of dissolved O, and 
O,/Ar ratio measurements in groundwater from monitoring wells were 
filled directly into 20 mL headspace vials allowing many volumes of 
overflow to displace air, poisoned with 100 uL of a saturated zinc 
chloride solution, quickly crimp-sealed with butyl septa (with no 
headspace bubble). At the laboratory, a 5mL headspace was estab- 
lished by replacement with ultra-high purity helium. Vials were equi- 
librated with headspace on a shaker table at room temperature for 
several days. Aliquots of the headspace (containing O3 and Ar from the 
water sample) were directly injected into the sample inlet system 
(including a 2m, 5A molecular sieve GC column) and routed into the 
isotope ratio mass spectrometer (Isoprime 100, multi-collector). 
Injections of lab air and lab air-equilibrated water were used as working 
standards. Based on repeat analyses the precision of the 680 mea- 
surements was +0.1%o and of the O2/Ar was +0.05. 


Cell staining, enumeration, and biomass yield estimations 

To fix cells we added 4 mL formaldehyde solution (37%) to a 50 mL 
groundwater sample (f. c. ~2.7%). The sample was stored at 4 °C until 
further processing. 1mL sample was diluted in 10 mL 1x phosphate 
buffered saline (PBS) and filtered onto a polycarbonate filter (0.1 um 
pore size, 25mm diameter, Millipore Sigma) using a mixed cellulose 
ester membrane filter (0.45 um pore size, Sartorius) as support. The 
filter was rinsed with 10 mL of 1x PBS, dried, and stored until further 
processing. A section of each filter was stained with 1ug mL DAPI 
(4’,6’-diamidino-2-phenylindole) for 10min at room temperature, 
washed with deionized water and 80% ethanol, and dried. Filter pieces 
were mounted on microscope slides using 4:1 Citifluor:Vectashield 
solution (VWR and Vector Laboratories) and stored at -20 °C. Cells 
were visualized and counted using an Axio Imager A2 (Zeiss, Jena, 
Germany) equipped with an X-Cite 120 LED (Excelitas, Waltham, USA) 
fluorescence light source and 12.5 x 12.5 mm ocular grid. -40 grids per 
filter were counted to obtain a robust dataset. To estimate the cell 
number that can be sustained per mol oxygen (Supplementary Data 8) 
we used published values for microbial cell content’’ and biomass 
yield”. 


16S rRNA gene library preparation and amplicon sequencing 
100-400 mL groundwater sample was centrifuged at 4000 x g for 1h 
at 4°C. The supernatant was discarded, the pellets combined in a 2 mL 


tube and stored at -80°C until further processing. (Note: We test 
collected cells of 8 wells also via filtration on 0.1 um pore size GTTP 
filters (Millipore Sigma), sequenced DNA, and compared the commu- 
nity structures to those derived of the pellets. The communities were 
nearly identical). Genomic DNA was extracted using the DNeasy Pow- 
erLyzer PowerSoil Kit (12855-100, QIAGEN) according to manu- 
facturer’s protocol with a minor modification; cells were lysed by bead 
beating at 4m s” for 45s using a Bead Ruptor 24 (OMNI). Extraction 
blanks were processed to detect potential laboratory contamination 
during extraction. DNA concentrations were measured fluorome- 
trically using a Qubit 2.0 (Thermo Fisher Scientific, Canada). The 
bacterial 16S rRNA gene v3-4 region was amplified using S-D-Bact-0341- 
a-S-17 (5’--CCTACGGGAGGCAGCAG-3’) and S-D-Bact-0785-a-A-19 (5- 
GACTACHVGGGTATCTAATCC-3‘)’®. The archaeal 16S rRNA gene v6-9 
region was amplified using the primer pair S-D-Arch-0915-a-S-20/S- 
*-Univ-1392-a-A-15 (5’-AGGAATTGGCGGGGGAGCAC-3’, 5’-ACGGGCGG 
TGTGTRC-3’)’*. PCRs consisted of 8 uL (1-10 ng) DNA template, 2.5 uL 
of each primer (f.c. 1M), 12.5 uL 2x Kapa HiFi HotStart Ready Mix 
(Kapa Biosystems, Wilmington, MA, USA) and PCR-grade water ad 
25 uL. For bacteria, a touchdown PCR program was used for improved 
annealing: initial denaturation at 95 °C for 3 min, 10 cycles of 95 °C for 
30 sec, 60 °C for 45 sec (touchdown -1°C per cycle), 72 °C for 60 sec, 
followed by 20 cycles of 95°C for 30sec, 55°C for 45sec, 72°C for 
60sec, and a final extension at 72°C for Smin. For archaea, the 
touchdown PCR program was: initial denaturation at 95 °C for 5 min, 10 
cycles of 95 °C for 30 sec, 62 °C for 45 sec (touchdown -1°C per cycle), 
72°C for 60sec, followed by 20 cycles of 95°C for 30 sec, 60°C for 
45 sec, 72 °C for 60 sec, and a final extension at 72 °C for 5 min. PCRs 
were performed in triplicate, pooled, and purified using 56 uL Agen- 
court AMPure XP beads (Beckman Coulter, Indianapolis, USA) per 
pooled PCR product (~65-75 uL) following manufacturer’s instruc- 
tions. Amplicons were indexed using 5 uL purified PCR product, 5 uM 
of each Index Primer (f.c. 1 uM), 25 uL 2x Kapa HiFi HotStart Ready Mix 
and 10 uL PCR-grade water. Indexing PCR program: initial denaturation 
at 95 °C for 3 min, 10 cycles of 95 °C for 30 sec, 55°C for 45 sec, 72°C 
for 60 sec, and a final extension at 72 °C for 5 min. Indexed amplicons 
were purified with Agencourt AMPure XP beads. The concentration 
and size of indexed amplicons were checked with a Qubit 2.0 fluo- 
rometer and Agilent 2100 Bioanalyzer system (Agilent Technologies, 
Mississauga, ON, Canada), respectively. Indexed amplicons were 
pooled in equimolar amounts and sequenced using IIlumina’s v3 600- 
cycle (paired-end, MS-102-3003) reagent kit on an Illumina MiSeq 
benchtop sequencer (Illumina Inc., San Diego, CA, USA) after all DNA 
extraction blanks and PCR reagent blanks were confirmed for negative 
amplification. 


Community analyses 

Raw amplicon sequences were analyzed using DADA2 v1.16”. Briefly, 
forward, and reverse reads were quality-trimmed to 275bp and 
205 bp, respectively, and primer sequences (17 bp forward, 21bp 
reverse) were removed. Reads with more than two expected errors 
were discarded, paired reads were merged, and chimeric sequences 
were removed. Species level taxonomy was assigned with the sil- 
va_nr_v138 train_set and silva species _assignment_v138. After quality 
control and the removal of blanks and technical replicates we 
obtained 64 archaeal and 110 bacterial amplicon datasets. Archaeal 
datasets contained a total of 3.89 x 10° sequence reads belonging to 
2633 unique ASVs. Archaeal samples had on average 6083 + 6387 
reads (mean + standard deviation) and 69 + 81 ASVs (Supplementary 
Data 10). Bacterial datasets comprised a total of 4.65 x 10° sequence 
reads belonging to 14665 unique ASVs. Bacterial samples had on 
average 4.23+1.57x10* reads and 272+221 unique ASVs (Supple- 
mentary Data 11). Supplementary Data 10 and 11 also include NCBI 
Project and Sample Accession numbers for each archaeal and bac- 
terial dataset. Archaeal and bacterial ASV nucleotide sequences and 
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taxonomy are listed in Supplementary Data 12 and 13, respectively. 
The archaeal or bacterial ASV-by-sample table (Supplementary 
Data 3, 4) was used to determine the number of observed ASVs, 
absolute singletons, relative singletons, relative abundance, and 
composition. Alpha diversity (richness, Shannon entropy, Inverse 
Simpson Diversity, and Chaol estimated richness) was calculated 
from the ASV-by-sample table using a subsampling approach to 
account for unequal sampling effort. We used 1008 and 6796 ran- 
domly chosen reads from each archaeal and bacterial sample, 
respectively. Differences in diversity between conditions were tested 
using the Wilcoxon signed rank-test (ggsignif) as implemented in 
geplot2°°. Bray-Curtis dissimilarities between all samples were cal- 
culated and used for two-dimensional nonmetric multidimensional 
scaling (NMDS) ordinations with 20 random starts. The further two 
samples (circles) are apart in the NMDS the more different are their 
underlying communities. Environmental parameters were tested 
with a partial Redundancy Analysis using Hellinger-transformed ASV 
data. All analyses were carried out with VisuaR v02 (https://github. 
com/EmilRuff/VisuaR) a workflow based on the R statistical envir- 
onment v4.1.0-v4.2.1, including among others the packages vegan*', 
ggplot2® as well as custom R scripts. The VisuaR workflow, the used 
packages, and versions are available in Supplementary Data 14. 


Metagenome sequencing and analysis 

Metagenomes were processed and sequenced at the Center for Health 
Genomics and Informatics in the Cumming School of Medicine, Uni- 
versity of Calgary. Genomic DNA was sheared into ~350 bp fragments 
using a S2 focused-ultrasonicator (Covaris, Woburn, MA). Libraries 
were prepared using the NEBNext Ultra II DNA Library Prep Kit for 
Illumina (E7103, New England Biolabs, Ipswich, MA) according to the 
manufacturer’s protocol, including size selection with SPRIselect 
magnetic beads (B23318, Beckman Coulter, Indianapolis, IN) and PCR 
enrichment (eight cycles) with NEBNext Multiplex Oligos for Illumina 
(E7335L, New England Biolabs, Ipswich, MA). DNA concentrations were 
estimated using qPCR and the Kapa Library Quantitation Assay for 
Illumina (Kapa Biosystems, Wilmington, MA). Genomic DNA was 
sequenced on an Illumina NovaSeq 600 sequencer (Illumina, San 
Diego, CA) using a 300 cycle (2 x 150 bp) S1 flow cell. Quality control 
was performed on raw, paired-end Illumina reads using BBDuk 
v38.90™, including trimming, filtration of contaminants, and clipping 
low-quality ends. Reads that passed quality control were assembled 
into contigs using Megahit v1.2.9% and contigs of <500 bp were not 
included in subsequent steps. Reads from each sample were mapped 
to each of the 5 samples using BBMap v38.90® and sequencing depth 
profiles were generated using the ‘jgi summarize_bam_contig depths’ 
function from MetaBat2 v2.15**. Assembled contigs were binned into 
metagenome-assembled-genomes (MAGs) using MetaBat2, CONCOCT 
v1.1.0% and MaxBin2 v2.2.7°°. DAS-Tool v1.1.3% was used to integrate 
MAGs produced by the three binning tools. The MAG properties 
including completeness and contamination were assessed by CheckM 
v1.2.0°°. To estimate the relative abundance of MAGs we determined 
the percentage of reads that mapped to each bin and the percentage 
that mapped to unbinned contigs using checkm coverage and checkm 
profile and then calculated (% reads mapped to bin) * (100 — (% reads 
mapped to unbinned contigs)). MAGs were classified using GTDB-tk 
(version 2.1.0, database release r207)°’. Metagenomic short reads were 
mapped to the SILVA SSU reference database v138”° to assign nearest 
taxonomic units, as well as full-length 16S/18S rRNA gene sequences 
were reconstructed from metagenomes using phyloFlash v3.4”! and 
compared to ASVs using blastn v2.12.0°°. Transfer RNA, ribosomal 
RNA, CRISPR elements, and protein-coding genes including nitric 
oxide dismutase (nod) and chlorite dismutase (c/d) coding genes were 
predicted and annotated using MetaErg v2.2.x”*. Amino acid sequen- 
ces on contigs were additionally searched against representative nod 


and cld sequences from NCBI using blastp v2.12.0 (e-value le 7°) and 
aligned using Clustal Omega v1.2.4". 


Contamination control 

To avoid and monitor contamination we have included controls at 
each step of the experiments. We have sterilized the sampling gear 
after each well, we included duplicate blind controls of regular water 
samples by using encrypted sample names, and we included blank 
controls for cell counting and sequencing. We have no indication that 
sample handling or DNA extraction introduced organisms that are 
known contaminants of low biomass samples” or DNA extraction 
kits’’. The storage in Nalgene bottles which could have had a potential 
enrichment effect, e.g., promoting growth of aerobic organisms after 
sampling, had no effect on community structure (Supplementary 
Fig. 11), as alpha diversity, beta diversity and composition did not show 
any differences when grouped into categories reflecting storage 
duration. 


Reporting summary 
Further information on research design is available in the Nature 
Portfolio Reporting Summary linked to this article. 


Data availability 

The archaeal and bacterial 16S rRNA amplicon data generated in this 
study have been deposited in the NCBI SRA archive under BioProject 
accession number PRJNA861683. The shotgun metagenomic data and 
metagenome-assembled genomes (MAGs) have been deposited in the 
SRA archive under BioProject accession number PRJNA700657. The 
comprehensive environmental data generated in this study have been 
deposited in the PANGAEA archive under accession number 952473”. 


Code availability 

The most recent version of the bioinformatic workflow VisuaR v02 
used to perform 16S rRNA gene-based community analyses is publicly 
available at Github (https://github.com/EmilRuff/VisuaR), the exact 
VisuaR analyses featured in this paper are deposited in Supplementary 
Data 14. 
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